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ABSTRACT 

The method originally developed by Kalnajs for the numerical linear stability analy- 
sis of round galactic disks is implemented in the regimes of non-analytic transformations 
between position space and angle-action space, and of vanishing growth rates. This al- 
lows effectively any physically plausible disk to be studied, rather than only those having 
analytic transformations into angle-action space which have formed the primary focus of 
attention to date. The transformations are constructed numerically using orbit integra- 
tions in real space, and the projections of orbit radial actions on a given potential density 
basis are Fourier transformed to obtain a dispersion relation in matrix form. Nyquist dia- 
grams are used to isolate modes growing faster than a given fraction of the typical orbital 
period and to assess how much extra mass would be required to reduce the growth rate 
of the fastest mode below this value. To verify the implementation, the fastest m = 2 
growth rates of the isochrone and the Kuzmin-Toomre disks are recovered, and the weaker 
m = 2 modes are computed. The evolution of those growth rates as a function of the halo 
mass is also calculated, and some m = 1 modes are derived as illustration. Algorithmic 
constraints on the method's scope are assessed and its application to observed disks is 
discussed. 
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1. Introduction 

Theoretical studies of the stability of thin stellar disks provide useful constraints on models of galactic 
formation and dynamics. Disk models can be excluded as unrealistic if they are found to be very unstable 
and the results of galaxy formation studies can be clarified where they lead to models close to marginal 
stability. When applied to observed disks, stability analysis provides a unique tool to probe what fraction of 
the mass is in luminous form via the requirement that there be enough extra matter to make the observed 
distribution stable over a period comparable to the age of the galaxy. 

Three main approaches have been adopted for disk stability analysis: direct N-body simulations, N-body 
simulations where the bodies are smeared onto a biorthonormal basis thereby solving Poisson's-equation im- 
plicitly, and linear modal analysis. The first has_been widely used, initially by Hohl (1971)13, and more 
recently, e.g. by Athanasoula & Sellwood (1985)tl. It provides a flexible tool of investigation which can 
be carried into the non-linear regime. These simulations, however, only provide insight into instabilities in 
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the statistical sense, and generally probe poorly the marginal stability regime. This drawback has recently 
been addressed by Earn & Sellwood (1995)0 who developed the second approach of solving Poisson's equa- 
tion through a biorthonormal basis. The stability of stellar systems has also been explored for spherical 
systems via a global "energy principle" (see Sygnet & Kandrup. (1984)Ej) but this appxqach has not been 
successfully implemented for disks because of resonances (Lynden-Bell & Ostriker (1967)113) . Moreover, the 
method usually provides stability statements which are of little practical use, since no time scale for the 
astrophysically relevant growth rates is available. _ _ _ 

The third method, linear modal analysis (Kalnajs (1977)^, Zang (1979)B Hunter (1992)liil) , is used 
here. In spite of the obvious limitation to perturbations of small amplitudes it has several potential ad- 
vantages. In this approach the integro-differential equation resulting from self-consistent solutions of the 
Boltzmann and Poisson equations is recast into a non-linear eigerbjaxthlerrj using an appropriately defined 
orthonormal basis. The method, pioneered by Kalnajs (1971-77)113113113 on the linear stability of galactic 
disks, involves the formulation of the dynamical perturbed equations in angle-action coordinates and the 
restriction to a (finiteL-set of perturbed densities which diagonalise Poisson's equation in position space. 
Kalnajs' original workli3 concerned the stability of the isochronc disk fat- which explicit transformations 
between angle-action and position-velocity coordinates exist. Zang (1979)123 studied the self-similar Mestel 
disks in the same way. More recently, the method was implemented by Hunter (1992)liil( via repeated eval- 
uation of elliptical integrals) for the infinite Kuzmin-Toomre disk. The approach of these authors involved 
calculating explicit transformations between angle-action variables and position-velocity variables and was 
therefore restricted to a very limited set of simple analytic .disks. More recently, .the relevant angle action 
integrals were computed by quadrature by Weinberg (1991)E3, and Bertin (1994)13 while studying the sta« 
bility of spherical systems for which the Hamiltonian is separable. Finally Vauterin and Dejonghe (1996)c3 
carried out a similar analysis for disks while performing the whole investigation in position space. 

Here, the first step of mapping the distribution function into angle-action space is implemented numer- 
ically by calculating the appropriate transformations from the results of integrating unperturbed orbits in 
the mean field of the galaxy. This then allows considerable freedom in the initial equilibrium to be studied, 
including the prospect of applying linear stability analysis to distribution functions recovered for observed 
galaxies. This method is not restricted to integrable potentials and could be generalised to 3D. It can in 
particular be implemented for measured distribution functions, where the potential is deduced from the 
rotation curve. 

In section 2 Kalnajs' matrix method is recovered following an indirect route corresponding to a rewriting 
of Boltzmann's and Poisson's equations directly in action space. The stability criteria are then reformulated 
in a form suitable for numerical evaluation of the matrix elements in Section 3. Details of the numerical 
method are presented and its range of validity is discussed. Section 4 contains results of calculations^for the 
isochrone and the Kuzmin-Toomre disks. These are in agreement with those found by Kalnajs (1978)t3, Earn 
& Sellwood (1995)0, Athanassoula & Sellwood (1985)0 and Hunter (1992)E1 Prospects and applications 
to observed data as probes of dark matter are sketched in the last section. 



2. Linear stability analysis revisited 

Linear stability analysis concerns the dynamical evolution of initially small perturbations in their own 
self-consistent field. The derivation (see eg Kalnajs (1977)ta, Sellwood & Wilkinson (1993)123) of the result- 
ing integro-differential equation is sketched here putting the emphasis on a coherent description in action 
space. This analysis is to be contrasted with that of Vauterin & Dejonghe (1996)c3 who chose to implement 
stability analysis for disks in position space alone. The resulting dispe rsion relati on, [Eq. (2.17)] , is fully equiv- 
alent to Kalnajs' Eq. K-15, but the intermediate integral equation, Eq. (2.lfj| , also—provides a connection 
with orbital stability analysis (Lynden-Bell (1979)0, Pichon & Lynden-Bell (1992)0 and Collett (1995)1) 
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2.1. The integral equation 

The Boltzmann-Vlasov equation, 



dF/dt + [H, F] = , 



(2.1) 



governs the dynamical evolution of an ensemble of collisionless stars. Here H is the Hamiltonian for the mo- 
tion of one star, F is the mass- weighted distribution function in phase space, and the square bracket denotes 
the Poisson bracket. Writing F = Fq + f, and W = ipp + ip, where Fo,ipo are the unperturbed distribution 



function and potential respectively, and linearising Eq. (2.1) in / and tp, the perturbed distribution function 
and potential, yields 

df/dt+[Ho,f]-bP,Fo}=0. (2.2) 

According to Jeans' theorem, the unperturbed equation, [Ho, Fq] = 0, is solved if Fq is any function only of 
the specific energy, s, and the specific angular momentum, h. 

Following Lynden-Bell and Kalnajs (1972)B, angle and action variables of the unperturbed Hamiltonian 
Hq are chosen here as canonical coordinates in phase-space. The unperturbed Hamilton equations are quite 
trivial in these variables, which makes them suitable for perturbation theory in order to study quasi-resonant 
orbits. The actions are defined in terms of the polar coordinates, (R, 9), by J = (Jr, Jg), where 



Jn 



(27T)- 



R dR , and Jg = h = R 2 



(2.3) 



Here [R] is a function of the radius, R, the specific energy, e, and the specific angular momentum, h, given by 
[R] = \j2e + 2ipo(R) — h 2 /R 2 . The angle between apocentres is = § h/R 2 [R]~ 1 dR. For a radial period 
Tr, and a pair of azimuthal and radial epicyclic frequencies fi, and k, given by 17 = Q/Tr and k = 2~k/Tr 
the phase-angles conjugate to Jr and h are tp = (tpR, (fig), where 



PR = « 







/" 


R 



dR , and 



(h/R 2 -n) R 



dR. 



(2.4) 



The stationary unperturbed Boltzmann equation Eq. (2.1) is solved by any distribution function of the 
form F = F (J), since [H , J] = 0. 

For growing instabilities, / and ip are both taken to be proportional to e~ , lu having a positive 
imaginary part. When expanded in Fourier series with respect to the angles tp, Eq. (2.2) becomes after 
simple algebra (Kalnajs 1977EJ) 

f m = m ^-^° /d \ m , (2.5) 



0/ — £l r 



where m is a integer vector with components (£, m), fig = m 1 m ■ fi = fi + In/m , fi stands for (k, fi), and 
fi p = uj/m. Here ip m and / m are the Fourier transforms of ip and / with respect to (<pn, <pg); for instance 



V> m (J) = J ip (R, 0) exp [—i (m • tp)] d 2 tp . 
Poisson's integral relates the potential, ip, to the density perturbation: 



ciR'.H') : 4;-G I |r R r/| dRdOdvRdvg 



(2.6) 



(2.7) 



This equation may be rewritten in order to make explicit the contribution from the interaction of orbits. Here 
again angle- action variables are useful, as a given unperturbed orbit is entirely specified by its actions. It is 
therefore straightforward to identify in Poisson's integral the contribution corresponding to the interaction 
of given orbits. Re-expressing this equation in terms of angles and actions (tp, J), using Parseval's theorem 
and taking its Fourier transform with respect to tp leads to: 



ip m (J) - 4ttG^ y / m , (J') A mm , (J, J') d 2 j' , 



(2.8) 
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where "0 m and / m are given by |Eq. (2.6)| and 

1 



(27T) 4 



exp^m' V - m • y) 2 , 2 
iR^BTi d¥ ' d¥ '- 



(2.9) 



The double sum in Eq. (2.8) extends in both £ and m from minus infinity to infinity where the radii R(<p, J) 
and R(<£>', J') are re-expressed as functions of these variables. Now |R — R'| depends on ip$ and Lp' ft in the 
combination ip' s — (pg = A(/> only. As |9 (93^ /<9 (A</5 <^e) | = 1, the order of integration in Eq. (2.9) may 
then be reversed, doing the ipg integratio n with Aip fixed. This yields 2irS mm i, so m' becomes (£', m) in the 

exp i (mAip ~ £'<p'ji + £^r) 



surviving terms. This gives for Eq. (2.9 



Ar, 



IR-R' 



-dA^d^ d(fin . 



(2.10) 



Each m mode therefore evolves independently. The dependence on m will be implicit from now on. So, for 
example, ipt is the Fourier transform of ip(R,9) with respect to both Lpn and ipg. 
Putting Eq. (2.5)| into Eq. (2.8) leads to the integral equation 



fa (Ji) 



J 



m 2 1 m 2 ■ dF Q /dJ 2 

Qi 2 ~ 



fa (J 2 ) d 2 J 2 



(2.11) 



This equation is the integral equation for the linear growing mode path an m-fold symmetry of a thin disk. 
Eq. (2.11)| was later approximated by Pichon & Lynden-Bell (1992)E3 Collett (1995)eJ Lynden-Bell (1995) to 
analyse the growth of the perturbation in terms of a Landau instability. 



2.2. The dispersion relation 



The perturbed distribution function and potential are now expanded over a potential-distribution basis 
{/ (n) }n, as 



/ (R, V)=Y, ««/ (n) (R, V) , and j> (R) = £ a n ^ (R) , 



(2.12) 



where the basis is assumed to satisfy P oisson's e quation Eq. (2.7), which when written in action space has 
Fourier components obeying (following Eq. (2.8) and given Eq. (2.10)): 



4 n) (3)=AnG^2 I /£° {3')Au> (J, J') d 2 J'. 

ci J 



(2.13) 



For basis functions scaling like exp(im'6) (preserving the axial symmetry described above), Eq. (2.6) applied 
to, say, ip( n \R) exp(im'6) gives for the Fourier mode £ : 



2tt 



(2.14) 



where Eq. (2.4) has been used to define the relative azimuthal increment 60((Pr) = tpe — 9. 

Expanding ip over this basis (using the same expansion for all £s) according to Eq. (2.12), inserting this 
expansion into Eq. (2.11) , multiplying by /j p ^*(J), integrating over J and summing over I yields 



m 2 1 m 2 -dF /dJ 2 An') 2 2 

*> nZ^K ^ (J 2 )dJ 2dJl 



(2.15) 
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Since fj?' belongs to the basis, it satisfies |Eq. (2.13)| and consequently 



E a » IE / # ^ ft ^ d2j = E ( E / 4T w ^ y 4:° d 2 J 2 



(2.16) 



Requiring that Eq. (2.16) has non trivial solutions in {a n } leads to the dispersion relation 

D{uj) = det |A - M (w) | = 0, (2.17) 



where the matrix M is defined in terms of its components (n',p) by the bracket on the r.h.s. of Eq. (2.16) 



while the matrix A corresponds to the identity if the basis (ip n . In) is bi-orthonormal, or to a matrix with 
components components (n,p) by the bracket pa the l.h.s. of Eq. (2.1 6)| otherwise. Equations |(2.17) was 
first derived in this context by Kalnajs (1977i1j) . In order to approximate the behaviour of a halo, a 
supplementary parameter q e]0, 1] is introduced 



D q (uj) = det | A - qM (w) | = , 



(2.18) 



as discussed below. 



2.3. Unstable modes & Nyquist diagrams 



The dispersion relations (2.1 



give, 
-iu)t- 



for each m, the criterion for the existence of exponentially growing 
■imO) . They are functions of a free complex parameter 10 corresponding 



unstable modes of the form exp( 
in its real part to m times the pattern speed of the growing mode, fi p , and in its imaginary part to the growth 
rate of the perturbation. The search for the growing modes is greatly facilitated by the use of Nyquist 
diagrams. Consider the complex to plane and a contour that traverses along the real axis and then closes 
around the circle at oo with Im(w) positive. The determinant D q is a continuous function of the complex 
variable to, so as lo traces out the closed contour, D q (oj) traces out a closed contour in the complex D q plane. 
If the D q (uj) contour encircles the origin, then there is a zero of the determinant inside the original contour 
and the system is unstable. In fact a simple argument of complex analysis shows that the number of loops 
around the origin corresponds to the number of poles above the line Im(w) = r\ — Constant. An intuitive 
picture of how the criterion operates on Eq. (2.17) is provided with the following thought experiment. 
Imagine turning up the strength of qG for the disk (i.e. turn on self-gravity, which is physically equivalent 
to allowing for the gravitational interaction to play its role; alternatively vary q, the ratio of relative halo 
support). Starting with small qG, all the M np are small, so for all lv's the determinant D q (u>) remains on a 
small contour close to det |A|. As qG is increased to its full value, either the D q contour passes through the 
origin to give a marginal instability or it does not. If it does not, then continuous change of qG does not 
modify the stability, and therefore the self-gravitating system has the same stability as in the zero qG case: 
it is stable. If, however, it crosses, and remains circling the origin then it has passed beyond the marginally 
stable case and is unstable. This picture is illustrated in the next sections where the image of the complex 
line mflp + irj, Q p g] — oo, oof, is plotted for various values of r\. 
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3. Numerical implementation 



3.1. Method 



Stability analysis, as implemented in the previous section, allows the treatment of a wide range of 
potentials and distribution functions which, in general, will not afford the explicit transformations between 
angle-action variables and velocity-position variables. 

Spline interpolation for tabulated potentials and distribution functions is used to compute the trans- 
formation of variables by numerically integrating orbits in the given potential. The whole calculation is 
performed on a grid of points in (Rq, Vq) space where Rq is the radius at apocentre, and Vq the correspond- 
ing tangential velocity. As may be seen from figure Fig. 3.4 in the next section, a 40 x 40 grid is adequate 
for the disks considered here. 

For each point in this grid, three orbits are calculated, one starting at the grid-point itself, and the 
others at small deviations in ho = RoVq, and along ho — constant. The angles and actions for these orbits 
are used later to calculate numerical derivatives of the distribution function. A fourth order Runge-Kutta 
integration scheme is employed for the orbits, stopping at the first pericentre for all but the most eccentric 
orbits where more than one oscillation is needed to give sufficient accuracy in the actions. Further points on 
the orbits corresponding to grid vertices are then calculated at exact subintervals of the orbital period. The 
Fourier components ipf" of a particular basis element along an orbit Eq. (2.14) are then given by: 



where t = <pr/k and 



4 n) (Ro,V ) = ±- e- 2 * m ' T « V (n) (R [*]) e lmSe ^dt , 
J fl Jo 



(3.1) 



59 ( m ) = k- 1 r ± dm - k- 1 ^ I A A VR = 9 (t) - (9 (t))t . 



(3.2) 



Numerically, this simply corresponds to taking a discrete Fourier transform (DFT) over the re-sampled 
points. The symmetry of an orbit means that although the argument is complex, the result is purely real. 
The matrix element defined in Eq. (2.16) is rewritten in terms of (Rq, Vq) and truncated in £: 



l{nU {dF/dh-m + dF/dJ ■(.) ,(»') 
\l ■ m + k ■ I — milp — ir] 



dJdh 



BRodVo 



di? dy 



(3.3) 



where the tp\ are now surfaces in (Rq, Vq) space given by the £'th order term of the DFT of the orbit at 

R , Vq. 

The inclusion of retrograde stars may be effected transparently by specifying a grid in Rq,Vq which 
includes negative Vq or, more efficiently, by noting that 



4 n) i-h) = ~ 

k -10 



cos (tip R (-h) - mS V (-hj) ^ (R [<p R i-h)])d<p R , 

= - [ cos {tipR (h) + mStp (h)) V (n) (R [(PR (h)})d<p R , 
* Jo 

= 4-l(h), 



(3.4a) 

(3.46) 
(3.4c) 



so with the appropriate sign switching the same set of orbits and Fourier modes may be used as for prograde 
stars. 

Having calculated the derivatives of the distribution function with respect to J and h by finite differences 
across the slightly displaced orbits described earlier and the Jacobian d(Jh)/d(Ro Vq) in the same manner, 
the problem of computing the matrix elements of equation Eq. (3.3)| reduces to one of integrating quotients 
of surfaces on the Rq, Vq grid. It is worth noting that, although it is the steps of the calculation up to 
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Figure 3.1: Nyquist diagram for the fastest growing mode of an isochrone/9 model 
for growth rates 77 = 0.12 (dashed curve), r\ = 0.14 (dot-dashed curve) and r\ = 0.16 
(dot-long-dashed curve). The magnification shows that the first two curves do enclose the 
origin whereas the third does not, indicating that the true growth rate is between 0.14 
and 0.16. 



this point which provide the generality of the implementation, their computational cost is relatively slight 
amounting in total, for example, to about a minute on a workstation with a specFp92 of 100. 

For vanishing growth-rates, the integrand in equation Eq. (3.3)| is a quotient of surfaces with a real 
numerator but a denominator which may vanish along a line within the region of integration. The integral 
then exists only as a principal value integral. Rather than attempting a brute force discretisation we approx- 
imate both numerator and denominator by continuous piecewise flat surfaces and employ analytic formulae, 
or their power series expansions, for each flat subsection. Each square of the Ro,Vq grid is treated as two 
triangles. According to the magnitudes and slopes of the two surfaces over a given triangle there are six dif- 
ferent approximations to be used for each component of the integral. For modes with finite growth-rates, the 
integrand in equation Eq. (3.3) is of a quotient of surfaces with a real numerator but a complex denominator 
and this prescription still holds. The extra generality afforded yields a solution which handles transparently 
the marginal stability case, where the resonances can be infinitely sharp. 
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Figure 3.2: Nyquist diagram for the fastest growing mode of an isochrone/9 disk with 
curves labelled as on 



Fig. 3.1. The dotted lines join points of the same ui. 
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Figure 3.3: Nyquist diagram for the fastest growing mode of an isochrone/9 as in 
Fig. 3.1, Fig. 3.2[ A logarithmic scaling has been applied near the origin mapping r — > 
log 10 (l + \Q a r)/a allowing the full topology to be seen more easily. 



The integration and summation then yields, for each ftp and 77, an n x n matrix M. Marginal stability 
corresponds to the vanishing of det|J — M\. As described in section 2.3 it is convenient to use Nyquist 
diagrams to locate the critical points: matrices corresponding to a hundred or so values of w are calculated 
for lines of constant r\. Providing the sampling is good enough, the number of times the line det| I — M\ 
encircles the origin in the complex plane gives the number of zeroes above the initial line mQ p + irj. An 
example of nyquist diagrams for the isochrone/9 disk (section 4) is illustrated in Fig. 3.1 and shown more 
quantitatively in Fig. 3.2| and Fig. 3.3] . The 3-dots-dashed curve corresponds to 77 = 0.16 and does not 
encircle the origin, (or rather it loops once clockwise and once anticlockwise) whereas the other two curves 
(77 = 0.12 and i] = 0.14) do, indicating that the growth rate for the fastest growing bisymmetric instability 
of this disk is between 0.14 and 0.16. 

This method is easy to apply to verify calculations where the growth rate and pattern speed are already 
known, but becomes more time-consuming when treating disks with unknown properties. Many automatic 
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schemes can be envisaged for locating zeroes in the complex plane but since computing the matrices is 
the slowest part of the calculation, we have found it most efficient to manage the search by hand with an 
interactive tool for examining Nyquist diagrams for different values of q (section 2.3). With a little experience 
the zeroes can be efficiently located even from poorly sampled diagrams with spurious loops which would 
easily lead to confusion in automated procedures. 



3.2. Validation 

There are two distinct types of errors to be assessed in this analysis: numerical errors arising from the 
discretization of integrals and from machine rounding; and truncation errors relating to the use of small 
finite bases, limited Fourier expansions and the calculation of only the inner parts of infinite disks. 

While in principle constraining the first type of errors is purely routine, it is worth noting that both 
the bases considered here require thejise of arbitrary precision languages (such as be) for their evaluation. 
As noted by Earn & Sellwood (1995)0, Kalnajs' basis requires the use of about fifty figures precision to get 
the first thirty basis elements to five figures. The problem, however, lies only in evaluating the polynomial 
terms because the coefficients are large and induce substantial cancelations. Qian's basis (1992)E3 not only 
requires extended precision to evaluate a single element but it is also extremely close to being singular. 

Besides evaluation of the basis elements, the only other part requiring special care is in the numerical 
derivatives of the distribution function with respect to J and h where a tight compromise must be made 
between the desire for a small interval to give an accurate derivative, and the loss of significance in the 
differences between angles and actions calculated for very close orbits. Nevertheless, a straightforward 
Runge-Kutta scheme for the orbits remains adequate in all the cases studied here. 



The second class of errors - those arising from the truncation of equations Eq. (2.17) - are best assessed 



by simply recomputing the same quantity with more basis functions, more Fourier harmonics, or a larger 



disk. The results of these tests for all the truncated quantities can be seen in Fig. 3.4 and Fig. 3.5. All these 
results are for the isochrone/9 model first studied by Kalnajs (1976)tia with the n — 7 series of Kalnajs' 
basis functions, which is further discussed in the next section. Since the calculation of the matrix M(w) in 



Eq. (2.18) is the most computationally intensive operation, the self-gravity parameter q is taken instead of 
the growth rate 77 as a convergence criterion. The q plotted in the figures is that which satisfies Eq. (2.18)| . 



It is found by bisection with respect to the winding number of the Nyquist diagram. A value 9=1 for 
the expected exact growth rate (given by Kalnajs (1978)Ej) is the asymptote towards which the calculation 
should converge when the appropriate range of parameters have been found. 

Fig. 3.4 shows how the sampling affects the result when the initial Rq,Vq grid has the same number 
of points in each direction. The different lines in this and subsequent figures indicate how q converges for 
different sizes of basis. The figure shows that for this disk, a sampling with a 40 x 40 grid reduces the 
sampling errors well below those from the basis truncation. The right-hand figure shows how the truncation 
of the disk influences convergence for different sizes of basis. If the instability is localised to the inner parts 
of the disk, including more of the outer regions should not affect the result. The behaviour seen here is due 
to stretching of the basis when i? ma x is increased, so more functions are required to sample the inner regions 
to the same resolution. If, however, enough functions are used, no change is seen beyond about i? ma x = 4.5. 

As remarked by Earn and Sellwood (1995)1 choosing an appropriate basis for a given instability has 
a signifi cant effect on the number of basis functions necessary to resolve the mode. This may be seen in 



Fig. 3.5 for Kalnajs' set of biorthonormal bases with 4 < k < 10. Ten functions with the k = 4 basis are 
required to give the same error in q as achieved with only six functions in the k = 10 basis. Kalnajs' remarks 
that some economies may be made by not using the same number of positive and negative radial Fourier 
modes but by centering the range on to, the order of the instability. This is illustrated in the right-hand 
figure, where the abscissa gives the number of radial modes on each side of a central value as labelled. The 
computational load is almost directly proportional to the number of radial modes to be used so this effect 
is well worth exploiting. Developing specific basis for different thin disks also speeds up the calculation, by 
reducing the number of functions required for convergence but it is not necessary. Any reasonable basis will 
get there in the end and the feasibility of the resulting computations considerably enhances the generality 
of this approach. 
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Figure 3.4: The self-gravity parameter, q as a function of the grid sampling, and the 
disk truncation. Data are shown for bases of various sizes as labelled. On the right, q 
is shown as a function of the truncation radius of the disk for different sizes of basis as 
labelled. More functions are needed for convergence with larger R because the inner, 
dominant, part of the disk is less well sampled when the basis is stretched over a larger 
domain. As with all the figures in this section, the disk is Kalnajs' isochrone/9 model. 
The basis is Kalnajs' biorthonormal set with k = 7. 
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Figure 3.5: The dependence of the self-gravity parameter, q, upon the order of the basis, 
the number of functions used in the calculation, and the number of radial harmonics. The 
left-hand figure shows q against the basis size for different bases from Kalnajs' set, labelled 
by the corresponding k parameter. On the right, 12 functions were used, with the number 
of radial harmonics on the abscissa taken on either side of a central value as labelled 
on the lines. Centering the harmonics on 2 (when looking for bi-symmetric instabilities) 
shows appreciably better convergence than taking equal numbers of positive and negative 
harmonics (centre 0). 
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Figure 3.6: The amplitude (left panel) and the phase (right panel) of components of the 
critical eigenvectors as a function of the number of basis functions. The rapid convergence 

Of the Shap^p^h&^d^&ppears clearly. Distribution Function 




Figure 4.1: left panel: isocontours of the mp = 4 distribution function in action space. 
The sharp break at zero momentum is an artifact of Kalnajs's trick which was used here 
to incorporate counter rotating stars; right panel: same isocontours for the prograde stars 
only in (Ra, Vo) - apocentre, velocity at apocentre - space. These variables are those used 
throughout the calculation to label the orbits. Note the envelope corresponding to the 
rotation curve of that galaxy. It appears clearly in this parametrisation that in this galaxy 
most orbits are of low ellipticity. 

4. Application: The isochrone and Kuzmin-Toomre disks 

The versatility of the algorithm for linear stability described above is now illustrated while recovering growth 
rates and pattern speeds of previously studied disk with old and new distribution functions, finding the modal 
response - in position space and in action space - of these disks both for bi-symmetric and lopsided modes. 
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Table 4.1: m = 2 growth rates and pattern speed of the isochrone/m^ model. 



The isochrone/mif model: bi-symmetric mode 


Model 




V 


6 


0.34 


0.075 


7 


0.37 


0.085 


8 


0.43 


0.125 


9 


0.47 


0.145 


10 


0.50 


0.170 


11 


0.53 


0.195 


12 


0.59 


0.210 



4.1. The Equilibria 

Three families of disks are studied here: two corresponding to equilibria for the isochrone disk (1961)0, 
and one for the Kuzmin-Toomre (1956)113 mass model. The isochrone disk is defined by its potential, -0 = 
GM/(b + rb) where rb 2 =R 2 + b 2 . The corresponding surface density is S = Mb{hi[(R + ri>)/b] — R/rb} 
/(2irR 3 ). The Kuzmin-Toomre potential is defined by ip = GM/rb .and the corresponding surface density 
is £ = M6/27rrf. Miyamoto (19X4)113, followed by Kalnajs (1976)13, Athanasoula & Sellwood (1985)13, 
and Pichon & Lynden Bell (1996)Eil chose specific forms of distribution functions assuming simple algebraic 
Ansatz for its expression in (e, h) space. Historically, these models were first used to construct families of 
maximally rotating disks {e.g. the Miyamoto/mM disk modelsU) while for othcr.s_ counter-rotating stars 
were re-introduced by simple, though somewhat arbitrary, tricks (Kalnajs (1978)113 described by Earn & 
Sellwood(19951I) . In this paper, The Kalnajs isochrone/m^ models are implemented in order to recover 
Kalnajs (1978)113 first linear stability results on differentially rotating disks. To demonstrate the flexibility 
of the method the truncated version of Hunter's Toomre-Kuzmin/mM models are also constructed, and the 
stability of the equilibria given by Pichon & Lynden Bell (1996) are analysed. 

Kalnajs's distribution family reads (in units of b and G = 1) 



where 



and 



g(x)=x 



f K (e, h) = 2 2m - 1 (V^2lh) m TP" 1 (-e) 

d(x m T m (x)) 



dx 



+ / (tx) m r m (tx) P^_ x (t) dt - |(m - l)mx m T m (x) 
o 



(x) 



log (r + VT 



VvT 



27rr 3 (_! _ v /TTr f )" 



(4.1) 
(4.2) 

(4.3) 



r=2x/(x 2 -l) 



Here P m stands here for the Legendre polynomial of order to. Pichon and Lynden-Bell's distribution family 
for the isochrone disk is given by 



fp (e, h) 



-1/2 



'-( — 

4tt 2 (m+1/2)!! ^2 \ds 



m+2 



{s 2 -l) m L{ S ) 



(4.4) 



s=l+fe 2 /2 



with L(s) = log(V s 2 — 1 + s) + Vs 2 — 1/s. Its contours for rap — 4 are illustrated in Fig. 4.1 while the Q 
profiles are given by Pichon & Lynden Bell (1996)cil( Fig PLB 5). Finally Miyamoto's distribution is 



, , . \ (2ro + 3) 2+2m / l -fts 

f M {e,h) = — — - — (-£) 2 Fi -to,-2-2to, -,— 
Z 7T \ lie 



(4.5) 



where 2-F1 is a terminating Hypergeometric function of the second kind. 
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Figure 4.2: the variation of the m = 2 growth rate (top left panel) and the pattern speed 
(top right panel) of an isochrone/mx family as a function of the self gravity parameter q 
and the Toomrc number Q. The bottom panels show the same results for the isochrone/mp 
family. The precision in the growth rates - at fixed truncation in the basis - drops for 
lower values of Q since these modes are more centrally concentrated. 



4.2. Characteristics of the linear wave 



4.2.1. Bi-symmetric m — 2 modes 

The pitch angle of the spiral response is commonly defined as cot(i) = (89 / d \og(R)) e where 9 = 9(H) 
at the crest of the spiral wave. In practice, it is best to calculate tan(i) = (dlog(R)/d9) R since R = R{9) is a 
bijection in [0, 2n/m\. Another useful set of quantities are defined by the radius, R1/2, (resp. the angle #1/2) at 
which the spiral response has decreased to half it maximum amplitude. The winding number ni/2 = #1/2/ (27r) 
yields a measure of the winding of the wave: the larger 711/2 the more wound it is. Comparing the position 
of the resonances and i? max to -R1/2 provides means to assess whether the truncation of the disk is likely to 
have generated spurious cavity waves. Specifically, it is required that R1/2 < Rolr < -Rmax so that the wave 
is well damped by the outer Lindblad resonance before the disk is truncated. Numerical simulation have 



suggested that R1/2 ~ -Rcor- Both statements are verified here for instance in Fig. 4.4 (resp. Fig. 4.5) which 




Figure 4.3: fits to the evolution of the m = 2 growth rate and the pattern speed 
of an isochrone/m/f (marked Kalnajs) and an isochrone/mp (makcd DLB) family as a 
function of a linear combination of the self gravity parameter q and the Toomre number 
Q. The left panel minimises the dispersion in r\ whereas the right panel that in w. The 
best fits are obtained for different linear combinations of Q and q. The solid line for 77 is 
r) = -0.05 - 0.25(Q - 2q) and that for w: w = 0.56 - 0.46(Q - q). 

Table 4.2: m = 2 growth rates and pattern speed of the isochrone/mp model. 



The isochrone/mp model: bi-symmetric mode 


Model 




V 


3 


0.21 


0.003 


4 


0.29 


0.032 


5 


0.35 


0.072 


6 


0.40 


0.105 



Table 4.3: m = 2 growth rates and pattern speed of the Toomre/wiM model. 



The Toomre/mM model: bi-symmetric mode 


Model 




V 


2 


0.598 


0.204 


3 


0.714 


0.294 


4 


0.810 


0.371 


5 


0.916 


0.445 



gives the linear response of a mx — 7 and a mx — 11 isochrone disk (resp. mx — 5 Kuzmin disk) for its 
first growing mode. These perturbations display the usual bar shaped central response with a loosely wound 
spiral response further out. 



Fig. 4.6 gives the evolution of cot(i), i?j 



/2, "-1/2- 



fi p , T), as a function of the 

Toomre number Q of the Kalnajs /mx disks. Table 4.1 - Table 4.3| give the growth rates and pattern speeds 
of the I rax /mp and jmM families. Note that for the new /mp family the bar mode given in Table 4.2 again 
have pattern speed well above the maximum of f2 — n/2 and hence do not display inner Lindblad resonance. 
Note also that the more compact Kuz'min -Toomre potential has larger growth rates and pattern speed as 
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Figirh-s- d-4: th p- i ri'crisitv response correspondiiig > *te4Jl£Jaste8t**m = 2 growing mode 
of the isochrone/7 and isochrone/11 model. Note t hat the colder /ll model yields a 
more tightly wound spiral as shown quantitatively in Fig. 4.6. This is expected since in 



the locally marginally radially unstable regime, the disk response should asymptotically 
match that of unstable rings. Its spiral response is also more centrally concentrated than 
that of its hotter counterpart. The solid circle corresponds to R1/2, the radius at which 
the wave has damped by a factor of two; the dotted circle to Corotation resonance; the 
dashed circle to the outer Lindblad resonance and the outer circle to R ma x = 5. 




Figure 4.5: same as* ]F4g^4 4| fnx^arTOyamoto- Hunter/3 model. Note that the response 
is much more centrally concentrated since the outer circle is R max = 3. This is expected 
since the Kuzmin-Toomre potential is more compact than the isochronc. 



expected since the dynamical time is shorter and the self gravity enhanced for those disks. Table 4.4 gives 
the growth rates for the second and the third growth rate of the ttik disks. Note that the growth rate of the 
second fastest growing mode of a tuk = 12 model was found by this method to equal uj + ii] = 0.461 + 0.145i, 
roughly within the error bars given by Earn & Sellwood (1995)0. Note that the number of radial nodes 
for these slower modes increase with the rank in stability (as illustrated in Fig. 4.1C for a m = 1 mode 
discussed in the next subsection) which is consistent insofar as winding decreases self gravity. Note also the 
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1 1.1 1.2 1.3 1.4 

Q 

Figure 4.6: the evolution of rj ,f2 p , cot(i), R\/2 and n 1 / 2 , as a function of the Toomre 
number Q of the isochrone/m^ disks. Note that the hotter disk have less wound spiral 
response. 

Table 4.4: second and third fastest m = 2 growing modes of the isochrone/mif model. 



The isochrone/mif model: bi-symmetric mode 


Model 


second mode 


third mode 


7 


0.29(4) + 0.04(l)i 


0.22(0) + 0.008(8)i 


10 


0.39(7) + 0.10(5)i 


0.27(3) + 0.039(5)i 


11 


0.42(8) + 0.12(8)i 


0.34(7) + 0.091(5)i 


12 


0.46(1) + 0.14(5)i 


0.26(4) + 0.051(2)i 



small relative jecror between the growth rates recovered for the Myamoto/Hunter models and those given by 
Hunter (1993)1111 for the Toomre-Kuzmin disks. Small residual discrepancies are expected given that Hunter's 
disks are infinite whereas those studied here are truncated at R = 5. [Fig. 4.2 gives for the ttik and mp 
families the evolution of the first growing mode as a function of the Toomre Q number (mass averaged in the 
inner region R < 2) of these disks and "the mass in the halo" as parametrised by q. An new attempted fit 
of a combination of Q and q is also given in Fig. 4.3 for both distributions; the relative dispersion illustrates 
the well known fact that the stability does not depend only on a simple combination of th£se numbers. 
These curves seem to be in qualitative agreement with those given by Vauterin & DejongheBl for different 
disk models. Here too, the pattern speed at fixed growth rate is a decreasing function of the Q number 
(dfl p /dQ n < 0). The conjecture of Athanasoula & SellwoodEJ of an asymptotic value of Q w 2 for marginal 
stability of fully self gravitating disk seems also consistent with these results. 
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Figure 4.7: the orbital r esponse in action space corresponding to the faste st growing 
mode described in Fig. 4.4 . The amplitude of the left hand side of Eq. (2.5) is plotted 
for the I = — 1 (ILR) fastest growing mode of an isochrone/12 model (left panel), and an 
isochrone/6 model, (right panel). Superimposed are the isocontours of the corresponding 
resonance il — re/2. The dashed line corresponds to the isocontour of a tenth of the pattern 
of the wave. Note that the hotter disk forms a Grand Design structure involving more 
eccentric orbits. 

4.2.2. Lopsided m = 1 modes 

The algorithm presented in section 3 is in principle inadequate to address the growth of m = 1 perturbations 



which for pure 
Following Zang 



^elf gravitating disks are forbidden since the perturbation does not conserve momentum. 
J, it is assumed here that the disk is embedded in some sufficiently massive halo to compen- 
sate its infinitesimal centre of mass shift. Fig. 4.£ gives the m = 1 modal response of the isochrone/mio,n 
disks while the corresponding pattern speeds and growth rates are given in Table 4.5. These modes have 



smaller growth rates than their bi-symmetric counterparts and are practically more difficult to isolate because 
they are close to other modes which grow almost as fast; the nyquist diagrams show many large l oops, al l 
of which must be properly sampled to avoid spurious contours encircling the origin as illustrated in Fig. 4.ij . 
Indeed, when investigating weaker and weaker modes the amount of looping in the corresponding nyquist 
diagrams should increase since modes always occur in pairs (u> ± rj) so that when r\ is small one is always in 
a regime corresponding at least to two (a growing and a decaying) modes. 

The physical mechanism leading to the appearance of these modes needs to be clarified, but presents 
little practical interest when their growth rates do not correspond to the fastest growing mode. The detailed 
analysis of m = 1 modes for the isochrone disk is therefore delayed until models of distribution functions 
which display weaker m = 2 modes are designed. 
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-5 5 -1.0 -0.5 0.0 0.5 1.0 

Figure 4.8: Nyquist diagram for the second m = 1 modes of the an isochrone/12 disk. 
There are more large loops than are seen for the bi-symmetric modes, and also more 
structure at small scales: on the right, the origin is shown enlarged under the mapping 
r — > (1 + 10 8 r)/8. All the space within the innermost of the three large concentric loops 
corresponds to the tiny cusp at the centre of the left hand panel. 




isochrone/10 and isochrone/11 model. The lopsided mode, which has a smaller growth 
rate, has a much more centrally concentrated response than its bi-symmetric counterpart. 



4.3. Application to observed disks 

The procedure described in section 3 makes no assumption on the nature of the distribution function or 
the potential of the disk. In particular, it is well adapted to distribution functions recovered from measured 
disks, where the radial derivative of the potential follows from the HI rotation curve of the disk while the 
distribution function itself is inverted from line of sight velocity profiles. The technique was tested on tables 
generated from Eq. (4.5)| . The agreement between the "theoretical" and "measured" (i.e. derived from a 
discretized representation of the distribution) growth rates was found to be in better than one percent. One 
peculiarity in this context lies in the finite difference chosen in the numerical derivatives of the distribution 
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Figure 4~M-ft L_ the m ^-^x density response corresponding to the second faster growing 
mode (ftp + irj = 0.18 + i 0.075) of the isochrone/12 model. The number of radial nodes in 
the response has increased by one compared to the corresponding faster mode illustrated 
on Fig. 4.9| for an isochrone/11 disk. 



Table 4.5: first and second growth rates and pattern speed of the m — 1 isochrone/mif 
model. 



The isochrone/mjc model: lopsided mode 


Model 


pattern speed 


growth rate 


9 


0.135 


0.0325 


10 


0.17 


0.065 


11 


0.20 


0.095 


12 


0.23 


0.125 



function with respect to J and h. While applying this analysis to observed data, special care should be taken 
in handling measurements relative to the core of the galaxy, since the relative fraction of counter rotating 
stars plays a crucial role in determining the growth rates of the disk, as pointed out already by Kalnajs 



and illustrated in Fig. 4.11. This appears clearly when recalling that counter-rotating stars add up to an 
effective azimuthal pressure in the inner core which therefore prevents the self gravity of the disk to build up. 
by orbit alignment. A non parametric inversion technique has been devised by Pichon & Thiebaut (1997)c3 
to recover the best distribution accounting for all the measured kinematics while handling specifically the 
counter-rotating stars. 
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Figure 4.11: the effect on q of artificially modifying the density of retrograde stars 
in Kalnajs' isochrone/9 distribution function. The retrograde part of the distribution 
function has been multiplied by the factor on the ordinate, keeping the prograde stars 
the same. Although in Kalnajs' model stars on retrograde orbits only account for a few 
percent of the total mass, they play a significant role in stabilising the disk. 

5. Conclusion & prospects 

A numerical investigation of linear stability of round galactic disks has proven successful in recovering 
known growth rates for the isochrone disk and the Kuzmin Toomre disks. The method is fast and versatile, 
and can be applied to realistic disks with arbitrary density and velocity profiles and relative "halo" support. 
Unstable bi-symmetric growing modes for new equilibria, and the fastest growing lopsided modes for Kalnajs' 
distribution function have also been isolated to demonstrate the code's versatility. The authors are currently 
applying the method in a study of the stability to families of disks parametrised by their relative temperature, 
compactness, fraction of counterrotating stars and halo support in order to identify the orbits responsible 
for the instability, and to probe the intrinsic (orbital or wave-like) nature of the m = 2 bar instability. The 
nature of lopsided m = 1 instabilities in disk galaxies is also under investigation. The implementation of 
linear stability analysis for observed galactic disk, yielding a direct relationship between the growth rate of 
the instability, the measured kinematical characteristics of the disk and the relative mass in the halo is also 
being investigated. The requirement of marginal stability will provide an estimate of the minimum halo 
mass. All ingredients will then be in place to probe the amount of dark matter required to stabilise observed 
galactic disks. 
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